                                                                                                         
% CALCULATE SOME MACRO VARIABLES
PD = sum(sum(PMAT.^(-theta).*Pdist));   % price dispersion
Cbar    = (wbar/chi)^(1/gam);               % consumption
Nbar    = Cbar*PD;                            % hours
U = Cbar^(1-gam)/(1-gam) - chi*Nbar;      % utility
Rss    = mu/beta;                                                    % monthly gross nominal interest rate
mbar   = nu*Cbar^gam/(1-1/Rss);  % nu*Cbar^gam/(1-1/Rss) in baseline %Cbar if mpy

Adjusters=lambda.*Pdist_sh;                                       % density of adjusting firms
vecAdjusters = Adjusters(:);

Pdistans = ones(nump,1)*pstar-Pgrid;                                 % matrix of distances from optimal price

AbsPdistans=abs(Pdistans);                                           % matrix of absolute price changes
Pincreases=Pdistans.*(Pdistans>eps^.5);                              % matrix of price increases
Pdecreases=Pdistans.*(Pdistans<-eps^.5);                             % matrix of price decreases

freqpchanges=sum(sum(Adjusters));                                    % percent of firms changing price in a period
massPriceIncreases=sum(sum(Adjusters.*(Pdistans>eps^.5)));           % mass of firms which increase their price   
massPriceDecreases=sum(sum(Adjusters.*(Pdistans<-eps^.5)));          % mass of firms which increase their price    

fracPriceDecr=1-massPriceIncreases/freqpchanges;                     % fraction of price decreases
fracPriceIncr = 1-fracPriceDecr;

EPchange=sum(sum(Pdistans.*Adjusters))/freqpchanges;                 % mean price change conditional on change
MeanAbsPchange=sum(sum(AbsPdistans.*Adjusters))/freqpchanges;        % mean abs price change conditional on change

EPincrease=sum(sum(Pincreases.*Adjusters))/massPriceIncreases;       % mean price increase
EPdecrease=sum(sum(Pdecreases.*Adjusters))/massPriceDecreases;       % mean price decrease
AvAbsDistFromPstar=sum(sum(abs(Pdistans).*Pdist));                   % average absolute distance from optimal price

StdPchanges=sqrt(sum(sum((Pdistans-EPchange).^2.*Adjusters))...
                               /freqpchanges);                       % standard deviation of price changes

KurtPchanges=(sum(sum((Pdistans-EPchange).^4.*Adjusters ))...
                               /freqpchanges)/StdPchanges^4;          % kurtosis of price changes

horizon = 24;     
hazd = hazard(horizon, R, Pdist, T2, P, lambda);

%%%%% prep histogram plots (low inflation)
lbound = min(Pgrid);
hbound = max(Pgrid);
numbins = length(Pgrid);
edges=[-inf linspace(lbound,hbound,numbins-1) inf];                         % boundaries of price change intervals

%   linspace(lbound,hbound,numbins)';
nbins=length(edges)-1;                                               % number of price bins
denspchanges=zeros(nbins,1);                                                 % initialize mass in each bin with zero
for i=1:nbins                                                        % density in nbins price change intervals
    denspchanges(i)=sum(vecAdjusters((Pdistans>edges(i)...
                            & Pdistans<=edges(i+1))))/freqpchanges;  % count of firms in price change interval k 
end

Pchangegrid = linspace(min(Pgrid), max(Pgrid), length(denspchanges))';


AverD=sum(sum(D.*Pdist));                                            % mean loss of all firms 
StdD=sqrt(sum(sum((D-AverD).^2.*Pdist)));                            % standard deviation of the loss of all firms 

vecD=D(:);                                                           % vetorizes matrix 
vecV=Vf(:);                                                           % vetorizes matrix 
vecPdist=Pdist(:);                                                   % vetorizes matrix 
vecPdisteroded=Pdist_sh(:);                                       % vetorizes matrix
vecAdjusters=Adjusters(:);                                           % vetorizes matrix
vecAbsPdistans=AbsPdistans(:);                                       % vetorizes matrix 
vecPincreases=Pincreases(:);                                         % vetorizes matrix 
vecPdecreases=Pdecreases(:);                                         % vetorizes matrix 

AbsPch_dens=[vecAbsPdistans vecPdist];                               % matrix with AbsPdistans and their density
AbsPch_dens=sortrows(AbsPch_dens,1);                                 % sort by AbsPdistans ascending 
MedAbsPdistans=AbsPch_dens(find(cumsum(AbsPch_dens(:,2))>=.5,1),1);  % median absolute distance from optimal price

D_density=[vecD vecPdist];                                           % matrix with D's and their density
D_density=sortrows(D_density,1);                                     % sort by D's ascending 
MedD=D_density(find(cumsum(D_density(:,2))>=.5,1),1);                % median D

V_density=[vecV vecPdist];                                           % matrix with V and their density
V_density=sortrows(V_density,1);                                     % sort by V ascending 
MedV=V_density(find(cumsum(V_density(:,2))>=.5,1),1);                % median V

Pinc_density=[vecPincreases  vecAdjusters/...
              massPriceIncreases.*(vecPincreases>eps^5)];            % matrix of price increases and their density
Pinc_density=sortrows(Pinc_density,1);                               % sort by price increases ascending 
MedPincrease=Pinc_density(find(cumsum(Pinc_density(:,2))>=.5,1),1);  % median  price increase

Pdec_density=[vecPdecreases vecAdjusters/...
              massPriceDecreases.*(vecPdecreases<-eps^5)];           % matrix of price decreases and their density
Pdec_density=sortrows(Pdec_density,1);                               % sort by price decreases ascending 
MedPdecrease=Pdec_density(find(cumsum(Pdec_density(:,2))>=.5,1),1);  % median price decrease

EPrice=sum((PMAT(:)).*Pdist(:));                                     % mean price (levels) (NOT D-S AGGREGATE)

ElogPrice=sum(log(PMAT(:)).*Pdist(:));                               % mean log-level price (NOT D-S AGGREGATE)

VarPrice=sum((log(PMAT(:))-ElogPrice).^2.*Pdist(:));                 % std of prices (levels)


NewPriceDeviations=ones(nump,1)*pstar-ones(nump,1)*log(EPrice);   % deviations of new prices from mean
NewPriceIncreases=NewPriceDeviations.*(NewPriceDeviations > eps^.5); % deviations of new price increases from mean
STDpincrease=sqrt(sum(sum((NewPriceIncreases).^2.*...         
    Adjusters.*(Pdistans>eps^.5)))/massPriceIncreases);              % standard deviation of price increases

absP_density=[vecAbsPdistans vecAdjusters/freqpchanges];             % absolute price changes and their density
absP_density=sortrows(absP_density,1);                               % sort by price changes ascending 
MedAbsPchange=absP_density(find(cumsum(absP_density(:,2))>=.5,1),1); % median absolute price change

fracSmallChanges=sum(vecAdjusters((abs(Pdistans)<=0.05)))...
                              /freqpchanges;                         % fraction of small price changes
fracSmallChanges25=sum(vecAdjusters((abs(Pdistans)<=0.025)))...
                              /freqpchanges;                         % fraction of small price changes

AvRevenue = sum(sum(Cbar*PMAT.^(1-theta).*Pdist));                 % Average (and total) firms' revenue

if exist('menucost','var')                                              % Menu Cost flow is freqpchanges*menucost*wbar as share of:             
   McostinRev=100*freqpchanges*menucost*wbar/AvRevenue;                 %  the flow of revenue of all firms
   McostinWageBill=100*freqpchanges*menucost*wbar/(Nbar*wbar);          %  the total wage bill 
   McostinCbar=100*freqpchanges*menucost/Cbar;                          %  consumption (units?)
end

AdjCostinRev = NaN;
TimingCostinRev = NaN;
LossinRev = NaN;

varlambda=sum(sum((lambda-freqpchanges).^2.*Pdist_sh));           % Cross-sectional variance of lambda   
CalvoMenuMetric=varlambda/(freqpchanges*(1-freqpchanges));           % Jim's measure of state-dependence

intensive = sum(sum(Pdistans.*Pdist_sh ));                        % average DESIRED price change
extensive = sum(sum(lambda.*Pdist_sh ));                          % average frequency of adjustment
selection = extensive*(EPchange-intensive);                          % selection effect
                                                                     % hazard histogram (based on Pdist_sh!)
edges=linspace(-eps^.5,1,101);                                       % boundaries of hazard intervals
npbins=length(edges)-1;                                              % number of hazard bins
density_lambda=zeros(npbins,1);                                      % initialize mass in each hazard bin with zero
for i=1:npbins                                                       % build density in hazard intervals
    density_lambda(i)=sum(vecPdisteroded((lambda>edges(i)...
                                         & lambda<=edges(i+1))));    % count of firms in hazard interval k 
end
npbins=length(density_lambda);                                       % number of p bins
npstep=(1/(npbins-1));
if abs(sum(density_lambda)-1)>eps^.5
   %sum(density_lambda)
   warning('Density L does not sum to one') 
   zero = NaN; % tell solver to check another point
   return;
end

[out, OPTindex] = max(Vf); 

optflexprice=log(theta/(theta-1)*wbar);               % log of optimal flexible price
Pmat_ergodicDistFlexPrice = Pmatrix(optflexprice, Pgrid, pstep);
ergodicDistFlexPrice = Pmat_ergodicDistFlexPrice.*(ones(nump,nump)*(Pdist));

profits = Cbar.*PMAT.^(-theta).*(PMAT-wbar);
revenue = Cbar.*(PMAT.^(1-theta));

averprof = sum(sum(profits.*Pdist));
averoptprof = sum(sum(profits.*ergodicDistFlexPrice));
avoptrev = sum(sum(revenue.*ergodicDistFlexPrice));

loss1 = (averoptprof - averprof)/ averoptprof;
loss2 = (averoptprof - averprof)/ avoptrev;

[Mstar, pstar] = get_pstar(Vf,Pgrid,mu);
Pdistans = ones(nump,1)*pstar-log(PMAT) ;   %EQUIVALENT but using new notation.

fracadjusters = sum(sum(lambda.*Pdist_sh ))   ;
desiredadj = sum(sum(Pdistans.*Pdist_sh )) ;
selection = sum(sum(Pdistans.*(lambda - fracadjusters).*Pdist_sh));
desiredinflation = sum(sum(Pdistans.*lambda.*Pdist_sh));

%Some additional moments
%Price-adjustment probability as a function of last period prices (assumes linear
%interpolation)

%Next period (log) prices as a function of last period prices (assumes linear
%interpolation)
Lambda = diag(lambda);
size_prime = R'*T2'*Lambda*(-Pdistans)./(R'*T2'*lambda);

interp_type = 'pchip';
dh = 1e-8;                              %window for integration
lambdaprime  = nan(length(Pdistans),24);
[intensive(1,1),grossExtensive(1,1),selection(1,1)] = calc_selection(Pdistans,Pdist_sh,lambda,D,model,interp_type,dh);
lambdaprime(:,1) = lambdaprimes(R'*T2',lambda,1);
for ii=2:25
    lambdaprime(:,ii) = lambdaprimes(R'*T2',lambda,ii);
    [intensive(1,ii),grossExtensive(1,ii),selection(1,ii)] = calc_selection(Pdistans,Pdist_sh,lambdaprime(:,ii),D,1,interp_type,dh);
end
%[intensive,grossExtensive,selection] = calc_selection(Pdistans,Pdist_sh,lambda,D,model,interp_type,dh);
sumAll = intensive+grossExtensive+selection;
relIntensive = intensive./sumAll;
relGrossExtensive = grossExtensive./sumAll;
relSelection = selection./sumAll;

lambda_prime = lambdaprime(:,1);  %R'*T2'*lambda;

